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Abstract  Laser  surgery,  or  laser-induced  thermal  ther- 
apy, is  a minimally  invasive  alternative  or  adjuvant  to 
surgical  resection  in  treating  tumors  embedded  in  vital 
organs  with  poorly  defined  boundaries.  Its  use,  however,  is 
limited  due  to  the  lack  of  precise  control  of  heating  and 
slow  rate  of  thermal  diffusion  in  the  tissue.  Nanoparticles, 
such  as  nanoshells,  can  act  as  intense  heat  absorbers  when 
they  are  injected  into  tumors.  These  nanoshells  can 
enhance  thermal  energy  deposition  into  target  regions  to 
improve  the  ability  for  destroying  larger  cancerous  tissue 
volumes  with  lower  thermal  doses.  The  goal  of  this  paper  is 
to  present  an  integrated  computer  model  using  a so-called 
nested-block  optimization  algorithm  to  simulate  laser  sur- 
gery and  provide  transient  temperature  field  predictions.  In 
particular,  this  algorithm  aims  to  capture  changes  in  optical 
and  thermal  properties  due  to  nanoshell  inclusion  and 
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tissue  property  variation  during  laser  surgery.  Numerical 
results  show  that  this  model  is  able  to  characterize  variation 
of  tissue  properties  for  laser  surgical  procedures  and 
predict  transient  temperature  fields  comparable  to  those 
measured  by  in  vivo  magnetic  resonance  temperature 
imaging  techniques.  Note  that  the  computational  approach 
presented  in  the  study  is  quite  general  and  can  be  applied  to 
other  types  of  nanoparticle  inclusions. 

Keywords  Laser-induced  thermal  therapy  • 
Nanoparticles  • Prostate  cancer  • Laser-tissue  interaction  • 
Bioheat  transfer  • Finite  element  method 


1 Introduction 

The  latest  statistics  shows  that  cancer  remains  one  of  the 
leading  causes  of  death  in  the  United  States  [22] . However, 
advances  in  nanotechnology  and  its  applications  in  bio- 
medical science  and  engineering  over  the  past  two  decades 
have  enabled  numerous  innovative  and  more  effective 
cancer  diagnosis  and  treatment  modalities  [3,  7,  10]. 
Treatment  by  traditional  surgical  resection  procedures  are  a 
surgeon’s  standard  approach  for  removal  of  well-defined 
primary  tumors  in  non- vital  regions.  This  technique  is 
extremely  invasive  and  usually  associated  with  high  mor- 
bidity. In  contrast,  thermal  therapies  employing  a variety 
of  heat  sources  including  laser,  focused  ultrasound,  and 
microwaves  have  benefits  over  conventional  cancer  treat- 
ment alternatives  [12,  13,  19,  25,  26].  These  treatment 
therapies  are  minimally  invasive  and  can  provide  an 
alternative  option  to  treat  solid  tumors  embedded  in  vital 
regions.  Technological  advancements,  such  as  actively 
cooled  applicators  and  high  power  diode  lasers,  have  made 
laser-induced  thermal  therapy  more  efficient,  economical, 
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and  safer  than  other  thermal  therapeutic  modalities.  Some 
advantages  include  that  laser-induced  thermal  therapy  can 
be  used  to  treat  tumors  more  rapidly  than  other  modalities 
and  have  more  control  over  perfusion  effects.  Additionally, 
lasers  do  not  require  a complicated  setup  that  involves 
grounding  pads  and  can  be  incorporated  safely  into  any 
imaging  environment,  including  MRI,  with  minimal 
induced  image  artifacts.  The  interaction  between  multiple 
probes  for  treating  larger  tumors  can  be  synergistic  and 
fully  compatible  with  MRI  for  online  monitoring. 

On  the  other  hand,  nanoshell-mediated  laser  surgery  is 
able  to  direct  thermal  energy  into  target  regions  delivered 
by  optical  fibers  to  provide  a lethal  dose  of  heat  while 
minimizing  damage  to  surrounding  tissue  [9].  In  particular, 
laser  surgery  promises  effective  treatment  of  small,  poorly- 
defined  metastases  or  other  tumors  embedded  within  vital 
regions.  In  this  study,  we  consider  a special  class  of 
nanoparticles  known  as  nanoshells,  which  can  act  as 
intense  infrared  absorbers  which  increase  the  thermal 
deposition  of  laser  energy  into  tumors.  In  particular, 
nanoshells  provide  a potential  means  to  (a)  enhance  the 
delivery  of  laser-induced  thermal  energy  via  distributing 
the  heat  source  from  the  fiber  to  the  surrounding  vascula- 
ture and/or,  (b)  provide  a highly  conformal  and  targeted 
approach  to  laser-induced  thermal  therapy  in  which  normal 
tissue  is  spared  and  tumor  tissue  is  ablated  with  a high  level 
of  specificity.  Typically  nanoshells  consist  of  a concentric 
spherical  dielectric  (silica)  core  and  a thin  metal  coating 
(Au)  shell.  The  diameters  of  nanoshells  are  usually  in  the 
110  to  120  nm  range  and  have  been  shown  to  be  effective 
mediated  agents  to  control  the  temperature  field.  Nano- 
shells possess  a highly  tunable  plasmon  resonance  which 
determines  the  particle’s  scattering  and  absorbing  proper- 
ties. The  plasmon  resonance,  one  of  the  nanoshell’s  optical 
properties,  can  be  tuned  across  a broad  range  of  the  light 
spectra  from  the  ultra-violet  to  the  infrared  by  controlling 
the  ratio  between  the  radius  of  the  core  and  the  thickness  of 
the  shell  layer  [1,  11].  When  nanoshells  are  injected  to  the 
target  region,  laser-induced  thermal  energy  can  be  deliv- 
ered to  specified  locations  and  greatly  enhance  heat 
absorption  in  tumor  regions  due  to  the  change  of  optical 
properties  in  the  tissue  [6]. 

To  design  optimal  nanoshell-mediated  laser  surgical 
protocols,  it  is  crucial  to  accurately  characterize  the  optical, 
thermal,  and  biological  response  of  tissues  to  therapies  [20, 
21].  The  major  challenge  is  that  these  properties  can  be 
difficult  to  measure  and  vary  over  time  during  the  treatment 
due  to  biological  alteration  in  tissues.  The  goal  of  this 
paper  is  to  present  a novel  nested-block  optimization 
algorithm  for  nonlinear  transient  bioheat  transfer  model  to 
simulate  laser  surgery  in  the  presence  of  nanoshells  and 
with  consideration  of  dynamic  changes  in  optical  and 
thermal  properties  of  the  tissue  due  to  biological  alteration. 


In  this  study,  a three-dimensional  finite  element  nonlinear 
transient  bioheat  transfer  model  with  input  of  laser-tissue 
interaction  calculation  from  Monte  Carlo  fluence  model  is 
constructed.  Numerical  results  show  that  this  model  can 
reliably  characterize  changes  in  tissue  properties  and 
accurately  predict  temperature  fields  comparable  to  those 
measured  by  in  vivo  magnetic  resonance  temperature 
imaging  (MRTI)  technique.  Although  the  validation 
experiments  are  conducted  for  treating  prostate  tumors 
inoculated  on  SCID  (severely  compromised  immuno-defi- 
cient)  mice,  the  computational  approach  presented  in  the 
study  is  quite  general  and  can  be  applied  to  other  types  of 
cancer  treatment.  Similar  optimization  strategies  are  also 
used  in  a dynamic  data-driven  framework  for  real-time 
surgical  control  [16]. 

2 Nanoshell-mediated  laser  surgery 

2.1  Tumor  preparation  and  inoculation 

Human  prostate  cancer  cells  (PC-3)  were  cultured  with 
HAM’s  F12  medium  with  10%  FBS  and  5%  penicillin- 
streptomycin.  Mice  were  anesthetized  with  isoflurane.  PC- 

3 cells  were  then  inoculated  in  the  backs  of  4-6  week  old 
SCID  mice  with  an  injection  volume  of  0.2  ml  of  5 x 106 
PC-3  cells  in  Matrigel  for  each  mouse.  Prostate  tumors 
were  grown  to  a tumor  burden  of  approximately  1.0  cc 
(equivalently  a spherical  tumor  with  radius  of  6.2  mm). 
Figure  la  depicts  a typical  PC-3  tumor  on  a SCID  mouse 
prepared  for  laser  irradiation. 

2.2  External  heating  and  nanoshell  inclusion 

Similar  to  previous  work  [9],  nanoshells  were  employed 
in  conjunction  with  extracorporeal  laser  irradiation  to 
enhance  thermal  deposition.  Nanoshells  were  injected  into 
the  mouse  tail  vein  24  h prior  to  laser  irradiation  to  enable 
adequate  accumulation  in  the  tumor  volume.  Nanoshells 
(furnished  by  Nanospectra  Biosciences  Inc.  Houston,  TX) 
are  composed  of  a silica  core  (with  a diameter  of  110  nm) 
and  an  outer  gold  shell  (thickness  of  15  nm).  The  shell 
geometry  was  optimized  to  enable  maximum  absorption 
at  wavelength  808  nm  to  enhance  thermal  deposition. 
Figure  lb  illustrates  the  distribution  and  geometry  of  the 
nanoshells  employed  in  laser  therapy  experiments. 

By  tuning  the  ratio  of  the  core  diameter  and  shell  thick- 
ness, nanoshells  can  be  made  to  absorb  or  scatter  light  at  a 
desired  wavelength  across  visible  and  near-infrared  (NIR) 
wavelengths.  This  optical  tunability  permits  optimal  design 
of  laser  surgical  protocols  with  a peak  optical  absorption  in 
the  NIR  region.  For  instance,  laser  surgery  for  deeper  tissue 
requires  light  in  the  NIR  region  where  tissue  has  the  highest 
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Fig.  1 a SCID  mouse  with  a prostate  tumor  (left),  b MRI  image  of  A-A  cross  section  with  illustrated  nanoshell  distribution  (right) 


transmissivity.  The  ability  of  gold-nanoshells  to  convert 
strongly  absorbed  light  into  localized  heat  can  be  exploited 
for  the  targeted  laser  therapy  of  cancer.  Thus,  effective  tar- 
geting of  nanoshell  bioconjugates  specifically  to  cancer  cells 
combined  with  the  high  absorption  cross-section  of  the 
nanoshells  in  the  laser  excitation  region,  generates  increased 
temperatures  sufficient  to  produce  irreversible  cell  and  tissue 
damage  to  subcutaneous  tumors  while  keeping  laser  energy 
at  a lower  level  so  that  cells  outside  of  the  target  region  are 
minimally  damaged. 

2.3  Temperature  measurement 

During  the  laser  surgical  experiments,  the  temperature  dis- 
tribution was  measured  by  in  vivo  MRTI  with  the  proton- 
resonance  frequency-shift  method  [8,  17].  All  experiments 
were  performed  at  the  University  of  Texas  M.D.  Anderson 
Cancer  Center  in  Houston,  Texas,  on  a 1.5-T  MR  scanner 
(Excite  HD,  GEHT,  Waukesha,  WI)  equipped  with  high- 
performance  gradients  (23  mT  m-1  maximum  amplitude 
and  120  T m-1  s-1  maximum  slew  rate)  and  fast  receiver 
hardware  (bandwidth,  ±500  MHz).  Mice  were  imaged  with 
a three-inch  receive-only  surface  coil  specially  designed  for 
small  animal  imaging.  Tr  and  T2-weighted  images  were 
used  to  plan  and  localize  the  treatment  by  verifying  the 
position  of  the  laser  fiber  relative  to  the  imaged  region  prior 
to  irradiation.  The  skin  over  the  tumor  was  swabbed  with 
PEG  diacrylate  (M*  600,  Sartomer,  West  Chester,  PA)  as  an 
index-matching  agent.  MRTI  was  performed  by  using  a 
complex  phase-difference  technique  with  a fast,  two- 
dimensional  RF-spoiled  gradient-recalled  echo  sequence 
(TR/TE  = 49.5/20  ms,  flip  angle  = 30°,  bandwidth  = 
9.62  kHz).  To  achieve  a five-second  per  image  scan  rate, 
a rectangular  field  of  view  (6x4  cm2)  with  256  x 86 
encoding  matrix  was  used.  Also  the  reduced  bandwidth  was 


employed  to  minimize  gradient  heating  limitations  at  this 
small  field  and  enhance  the  signal-to-noise  ratio.  The 
acquired  voxel  size  was  1.07  x 0.82  x 3 mm3.  The  change 
in  temperature  from  baseline  after  a number  of  images  was 
extrapolated  from  the  complex- valued  MRTI  data  by  using 
the  temperature  dependence  of  the  proton  resonance  fre- 
quency shift  and  a temperature  sensitivity  with  coefficient  of 
—0.0097  ppm/°C  [8].  The  temperature  resolution  for  MRTI 
measured  temperature  difference  is  less  than  1°C.  Figure  2 
shows  typical  MRI  and  MRTI  images  along  with  a quantified 
image  that  combines  the  two. 

3 Mathematical  and  computational  models 

Simulating  laser  surgery  and  making  reliable  predictions  of 
the  temperature  field  requires  two  major  modeling  com- 
ponents: a bioheat  transfer  model  for  the  tissue  and  a laser 
source  term  that  characterizes  thermal  energy  deposited 
into  the  tissue.  Since  the  optical  properties  of  nanoshells 
can  be  designed  by  adjusting  the  ratio  between  the  diameter 
of  the  silica  core  and  the  thickness  of  the  gold  shell,  tissue 
properties  such  as  the  absorption  and  scattering  coefficients 
can  be  controlled  then  with  inclusion  of  nanoshells.  In  this 
section,  we  discuss  the  mathematical  and  computational 
models  for  bioheat  transfer  and  laser-tissue  interaction. 

3.1  Bioheat  transfer  model 

The  mathematical  representation  of  the  temperature  dis- 
tribution in  the  tissue  incorporates  both  the  Pennes  bioheat 
equation  for  the  thermal  effects  of  local  blood  perfusion 
and  an  expression  for  laser  energy  as  a thermal  source.  We 
consider  the  case  of  external  heating  provided  by  a laser 
source  on  the  surface  of  the  tumor. 
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Fig.  2 a MRI  image  of  a mouse  tumor,  field  of  view  is  4 x 6 cm2  and  thickness  associated  with  the  image  is  3 mm.  b a thermal  image  of  the 
heating  of  the  mouse  tumor  overlaid  onto  the  MRI  image,  c a quantified  thermal  image  of  the  mouse  tumor  with  color  scale  from  308  to  334  K 


Let  Q be  a bounded  domain  (tumor  region)  in  7 Z3  with 
r = 0QC  |J  0Qn  denoting  a Lipschitz  continuous  bound- 
ary. Equation  (1)  is  based  on  Pennes  bioheat  equation  [18] 
with  temperature  dependent  coefficients: 

pcP  ^ - V • (k(T)VT)  + co(T)c„(T  - Ta)  = Q(x,  t)  (1) 

in  Q 

The  thermal  conductivity,  k [J  s-1  m-1  K_1],  and  blood 
perfusivity,  co  [kg  s1  m-3],  are  assumed  to  be  bounded 
functions  of  the  temperature  field,  T = T(x,  t),  where 

k(T)  = k0  + k\  atan(k2(T  — fc3)) 
co(T)  = coo  + co\  atan(co2 (T  — m3)) 

and  k0[J  s— 1m— 1 K~l],kx[ J s_1m_1K_1], k2[K_1], k3 [K], 
mo[kg  s_1  m-3],  mi[kg  s_1  m-3],  co2[K_1],  and  m3[^T | are 
parameters  of  the  diffusivity  and  perfusivity  coefficient 
functions  defined  above.  The  specific  heat  of  the  tissue  and 
blood  are  given  by  cp  and  cb  respectively,  Ta  is  the  arterial 
temperature,  and  p is  the  density  of  the  tissue.  On  the 
Cauchy  boundary, 

-k{T)VT  • n = h(T  - T^)  on  0QC 

The  coefficient  of  cooling  is  denoted  h.  T ^ denotes  the 
ambient  temperature.  Q is  the  prescribed  heat  flux  on  the 
Neumann  boundary. 

-k(T)VT  n = Q on  0QN 

The  temperature  field  is  propagated  forward  in  time  from  a 
given  initial  condition,  T0. 

r(x,  0)  = T0  in  Q 
3.2  Laser-heating  model 

Laser-tissue  interaction  is  a complex  phenomenon  and 
usually  divided  into  optical  and  thermal  responses  [15,  24]. 
Typically,  laser  light  is  absorbed  by  tissue  and  converted 


into  heat.  A commonly  used  model  that  characterizes  the 
absorbed  heat  can  be  represented  by  the  rate  of  heat  gen- 
eration Q(x , t)  that  is  defined  as 


2(x,  t)  = na  <D(x,  t)  = nantrP(t) 


3exp(  — jUeff  llx  — Xo|| 


47t||x  — xo|| 

where  n,r  = Pa  + ~ 8),  tks  = \/3/VV 


(2) 


and  0(x,  t)  is  the  fluence  that  defines  the  amount  of  energy, 
in  the  form  of  photons,  passing  through  a unit  area  at  a 
point  in  space  per  unit  time.  P(t)  is  the  laser  power  at  time 
t.  The  parameters  pa,  ps  are  absorption  and  scattering 
coefficients  that  relate  to  laser  wavelength.  They  can  be 
considered  as  the  probability  of  absorption  and  scattering, 
respectively,  of  photons  in  tissues.  The  constant  g is  the 
anisotropic  factor,  and  x0  is  the  position  of  the  laser  source. 

Derived  from  the  light  diffusion  approximation  for 
monoenergetic  neutral  particles  by  solving  the  transport 
equation  [23],  (2)  is  valid  when  the  scattering  coefficient  is 
much  larger  than  the  absorption  coefficient  (ps  pa), 
which  is  the  case  in  most  of  the  soft  tissues.  The  basic 
assumption  in  our  case  is  that  the  light  is  emitted  from  a 
single  point  in  a diffuse  and  isotropic  manner.  In  addition, 
this  model  only  accounts  for  scattered  light,  i.e.,  the  pri- 
mary light  is  ignored.  This  assumption  is  reasonable  if  the 
region  of  interest  is  far  from  the  source.  Lor  external 
heating,  these  assumptions  may  be  violated.  In  our  case, 
the  light  from  the  laser  beam  was  collimated  as  it  enters  the 
tissue  and  is  incident  against  the  skin  with  a flat  beam  spot 
between  0.5  mm  and  1.0  cm  in  diameter. 

In  this  study,  we  use  both  the  classical  isotropic 
approach,  (2),  and  the  Monte  Carlo  method  to  solve  for  the 
fluence  (see  e.g.  [23]).  Generally  speaking,  the  Monte 
Carlo  method  is  considered  to  be  more  accurate  but  com- 
putationally more  expensive  than  the  classical  approach. 
We  use  the  Monte  Carlo  method  here  to  verify  the  accuracy 
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of  the  classical  isotropic  approach.  The  two  major 
assumptions  that  are  used  in  the  Monte  Carlo  model  are:  (a) 
the  photons  will  be  distributed  in  a cylindrically  symmetric 
manner  with  the  initial  direction  of  the  beam  as  the  axial 
direction,  and  (b)  the  probability  distributions  being  used 
for  the  scattering  angles  and  mean  free  path  length  are 
known.  To  compute  the  fluence,  we  consider  the  following 
three  quantities  of  interest  as  usual. 

1.  the  absorption  coefficient,  \ia,  the  average  number  of 
photons  absorbed  per  unit  length, 

2.  the  scattering  coefficient,  fis , the  average  number  of 
photons  scattered  per  unit  length,  and 

3.  the  anisotropic  factor,  g,  the  expected  value  of  the 
cosine  of  the  deflection  angle. 


All  three  quantities  depend  on  light  wavelength,  tem- 
perature, and  space,  due  to  heterogeneity  of  the  tissue. 
However,  the  space  dependency  is  usually  ignored. 
Thus,  fia,  fis,  and  g are  typically  given  as  functions  of  light 
wavelength  and  temperature  only.  It  has  been  found  that 
values  of  these  three  quantities  remain  fairly  constant  for 
temperatures  below  70°C  [14]. 

The  Monte  Carlo  simulation  for  laser  irradiation  starts 
with  following  the  random  paths  of  many  individual  pho- 
tons by  sampling  probability  distributions.  The  main  idea  is 
described  as  follows.  The  algorithm  starts  by  initiating  a 
photon  with  a weight  of  W = 1 and  a position  and  direction 
in  space.  A random  number  between  0 and  1 is  then  gen- 
erated and  used  as  a value  of  a mean  free  path  length  by 
correlating  it  to  the  associated  probability  distribution.  The 
photon  is  then  moved  with  this  length  in  its  current 
direction.  Then,  the  algorithm  checks  if  the  photon  is  still 
in  the  tissue.  If  it  is  not,  the  internal  reflection  for  this 
particular  photon  is  done.  Otherwise,  a portion  of  the 
photon,  W • is  assumed  to  be  absorbed  at  that  loca- 

tion, the  weight  is  updated  and  then  two  more  random 
numbers  are  generated  to  obtain  the  new  azimuthal  angle 
and  the  deflection  angle  by  again  correlating  the  random 
numbers  to  the  respective  probability  distribution.  This 
process  is  repeated  until  the  photon’s  weight  has  reached 
some  cutoff  weight.  To  conserve  energy,  a roulette  scheme 
is  employed  here  to  decide  if  the  photon  should  be  dis- 
carded or  left  in  with  an  update  of  its  weight. 

The  output  of  Monte  Carlo  algorithm  is  the  number  of 
photons  absorbed  in  each  grid  cell  descretized  in  an  axi- 
symmetric  r— z plane.  Then,  the  fluence  can  be  obtained  by 


0>[z,  r] 


grid for] 

M ' V[z,  A ' Va 


where  grid[z,  r ] is  the  number  of  photons  absorbed  in  the  grid 
cell  [z,  r],  M is  the  total  number  of  photons  put  into  the 
simulation,  V[z,  r]  is  the  volume  represented  by  the  grid  cell 
[z,  r],  and  fia  is  the  absorption  coefficient.  The  grid  size  can 


be  refined  according  to  desired  tolerance.  Having  computed 
O [z,  r ] at  each  discrete  grid  cell,  it  will  be  projected  to  the 
three-dimensional  finite  element  mesh  and  used  to  calculate 
the  heat  generation  term  Q(x , t)  = ^a0(x,  t)  in  (1). 

3.3  Nested-block  optimization  algorithm 

In  this  section,  we  present  a nested-optimization  algorithm 
that  applies  to  the  Pennes  bioheat  transfer  model  [18].  The 
goal  is  to  capture  dynamic  changes  in  tissue  properties  due  to 
biological  alteration  as  well  as  nanoshell  inclusion.  In  order 
to  achieve  this  goal,  the  objective  function,  (5),  for  the  cal- 
ibration problem  is  defined  as  the  difference  in  space-time 
norm  between  the  computed  temperature  field,  T(x,  t),  and 
the  temperature  field  obtained  from  the  MRTI  experiments, 
Texp(x,  t),  integrated  over  the  entire  biological  domain  Q and 
time  duration  of  interest  [0,  t].  The  computed  temperature 
field  T(x,  t)  is  an  implicit  function  of  the  bioheat  transfer 
model  parameters  /?,  a vector  that  consists  of  the  location  of 
the  laser  probe,  absorption  and  scattering  coefficients,  and 
parameters  in  conductivity  and  perfusion  functions. 

The  main  question  to  be  addressed  in  this  section  is  how 
to  approximate  the  full  optimization  problem  efficiently  so 
that  the  prediction  can  be  made  prior  to  the  next  data 
arrival.  In  other  words,  the  speed  of  computation  needs  to 
be  faster  than  the  rate  of  data  acquisition,  at  least  faster 
than  the  time  duration  of  the  operation.  The  optimization 
problem  can  be  formally  stated  as: 

Find  such  that  ft)  is  minimized,  where 

T 

Q(m,p)=\J  J(m(x,t)-Texv(x,t))2<hdt  (3) 

0 Q 

In  order  to  speed  up  solution  time  for  the  approximation, 
(3)  is  replaced  by  a sequence  of  smaller  problems  with  less 
historical  data  in  the  time  dimension. 

Find  ft*i  such  that  Qi(T(P),P)  is  minimized,  where 

Tj  <T 

Qi{T{Pi),  ^ = \J  I ( m)(x , t)  - rexp(x,  t)f  dxdr 

0 Q 

(4) 

using  a smaller  time  window  [0,  tJ  C [0,  Ti+1]  C [0,  t], 
i = 1,2,3,....  Figure  3 illustrates  the  nested-block  optimi- 
zation process:  Calibration  — ► Prediction  — ► Validation/ 
Calibration  ->  Prediction,  as  in  vivo  MRTI  measurement 
data  continue  to  be  imported  into  the  model. 

Based  on  the  experiment,  the  temperature  field  is  cali- 
brated with  respect  to  MRTI  thermal  images  of  heating  of  a 
mouse  tumor  treated  with  nanoshells.  A nonlinear  form  of 
Pennes  bioheat  transfer  model,  ( 1),  is  calibrated  to  predict  the 
heat  transfer  in  tissues.  There  is  a trade-off  between  accuracy 
and  efficiency.  An  efficient  optimization  process  will  allow 
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Fig.  3 Illustration  of  nested-block  optimization  algorithm 


and  finite  difference  in  time  using  the  Crank-Nicolson 
scheme. 

The  laser  source  model  was  implemented  using  both  the 
isotropic  analytical  and  the  Monte  Carlo  models  as  described 
previously.  For  the  Monte  Carlo  calculation,  results  were 
obtained  by  averaging  over  10  runs  with  1 ,000,000  photons  in 
each  run.  Parameters  used  were  pa  = 2.15  cm-1  and  ps  = 
14.2  cm-1.  The  grid  size  was  chosen  A r = 0.015  cm-1  and 
Az  = 0.005  cm-1.  The  cumulative  distribution  function  for 
the  mean  free  path  length  is  defined  as 

P{s<Sl}  = (5) 


faster  solution  time  for  the  calibration  so  that  it  can  be  used  in 
real-time  surgical  control.  However,  time  constraints  of  real- 
time computation  do  not  permit  the  computation  of  the 
Hessian  of  the  objective  function,  which  would  allow  a more 
accurate  measure  of  rate  change  with  respect  to  each  model 
parameter  as  a solution  is  approached.  Thus,  we  compute  the 
gradient  of  the  objective  function  using  a limited-memory 
variable  metric  by  a quasi-Newton  optimization  method  [2]. 
As  shown  in  [16],  the  gradient  vector  of  the  quantity  of 
interest  can  be  written  as 


where  fit  = pa  + ps.  We  assume  that  the  probability 
distribution  of  the  azimuthal  angle  is  uniform.  The 
probability  distribution  of  the  cosine  of  the  deflection 
angle  is  defined  as 


/(cos  9)  = 


1 ~g2 

2(1  + g2  — 2gcos  9 )3^2 


(6) 


which  is  known  as  the  Heyney-Greenstein  scattering 
function,  where  the  anisotropic  factor  g is  set  to  0.9  (a 
typical  value  for  soft  tissues). 


VG(  P)  = 


[VT  • Vp],  [atan(fc2  ( T - k3))  VP  • Vp], 

, [(T  - Ta)  -p],  [atan(co2(7’  - co3))(r  - Ta)p], 


ki  ( T-k 3)  Vr-Vp 
1 +k\  (T-k3)2 


—k\  ki  VT-Vp 
1 +k\  ( T-k,)2 


0)1  ( T-CD3)(T-Ta)p 

-coi  co2(T—Ta)p 

1 +(D22(T-(D'$)2 

1+cd22{T-co3)2 

> 

l 


where  [ ☆ ] denotes  fQ  [☆]  dx  dt  and  p(x , t)  is  the  solution 
to  the  adjoint  problem 

-PcP^-V-mVp)+^(T,P)VT-Vp  + CD(T,P)p 
+ ^(T,p)P(T-Ta)  = T-Texp  in Q 

with  the  boundary  conditions 

—k(T)\7p  ■ n = h p on  0QC  — k(T)Vp  • n = 0 
on  6Q^v 

and  the  terminal  condition 
p(x,  t)  = 0 in  Q 

3.4  Computational  implementation 

The  major  computational  modules  and  data  flow  for  the 
integrated  nonlinear  transient  bioheat  transfer  computa- 
tional model  with  laser  heating  are  shown  in  Fig.  4. 

The  nonlinear  transient  bioheat  transfer  model  is 
implemented  using  a finite  element  discretization  in  space 


Both  bioheat  transfer  and  optimization  modules  are 
implemented  in  parallel  using  an  h—p  finite  element 
method  built  upon  an  in-house  code  and  can  be  executed 
efficiently  on  a multi-processor  machine.  The  framework 
of  this  adaptive  finite  element  code  has  been  developed  and 
thoroughly  tested  over  last  two  decades  [4].  Computations 
were  carried  out  using  Linux  clusters  hosted  at  the  Texas 
Advanced  Computing  Center  at  Austin.  Each  computing 
node  has  two  Xeon  Intel  Duo-core  64-bit/2.66GHz  pro- 
cessors. Using  up  to  50  processors,  we  are  able  to  make 
real-time  predictions  two  minutes  ahead,  which  only  takes 
10  s computing  time. 

4 Results 

Data  acquired  for  calibration  are  shown  in  Fig.  2.  Fig- 
ure 2a  is  a 49  x 56  pixel  MRI  image  of  a mouse  tumor. 
The  field  of  view  is  4 x 6 cm2  and  the  thickness  associated 
with  the  image  is  3 mm.  An  external  heat  source  was 
applied  to  a mouse  tumor.  Sixty  thermal  images  were 
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Fig.  4 Block  diagram  of  major 
computation  modules 


C.ili  bi.ition  Feed  tuck  Loop 


acquired  at  every  5 s.  A single  time  instance  of  the  thermal 
images  is  shown  in  Fig.  2c.  An  overlay  of  the  thermal 
image  onto  the  geometry  is  shown  in  Fig.  2b. 

In  order  to  compare  computational  results  with  the 
experimental  data,  the  in  vivo  temperature  measurement 
(MRTI)  data  are  interpolated  onto  a finite  element  mesh  of 
the  biological  domain  generated  from  the  MRI  data.  In 
Fig.  5,  the  results  at  time  t = 300  s are  shown.  The  initial 
mesh  consists  of  6,076  linear  elements  with  a total  of  8,000 
degrees  of  freedom.  The  MRTI  data  was  filtered  and  then 
projected  onto  the  finite  element  mesh,  which  is  shown  in 
Fig.  5a.  Figure  5b  compares  a cut-line  through  the  mid- 
plane of  the  unfiltered  thermal  image  to  the  filtered  thermal 
image.  The  Crank-Nicolson  scheme  was  used  to  propagate 
the  bioheat  transfer  equations  forward  in  time.  Figure  5c 
and  e show  the  thermal  contour  of  the  FEM  prediction 
corresponding  to  the  isotropic  laser  source  term,  (2),  and 
the  Monte  Carlo  input  respectively.  Parameters  used  in  the 
simulation  are  given  in  Tables  1 and  2.  The  initial  values  at 
the  beginning  of  the  calibration  process  are  taken  from  [5]. 

Figure  5d  and  f illustrate  a cut-line  comparison  of  the 
MRTI  data  to  the  FEM  prediction  using  the  isotropic  laser 
source  term  and  the  Monte  Carlo  source  term.  The  FEM 
model  with  uncalibrated  parameters  seems  to  produce 
inaccurate  results  compared  with  MRTI  data,  which 
strongly  indicates  that  the  tissue  properties  are  changed  due 
to  the  introduction  of  nano  shells. 

An  early  study  [6]  shows  that  the  model  of  bioheat 
transfer  with  a laser  heating  source  term  is  anticipated  to  be 
relatively  less  sensitive  to  the  attenuation  coefficients  of  the 
laser  source  term  as  compared  to  the  thermal  conductivity 
and  blood  perfusivity  coefficients.  Thus,  more  changes  in 
the  blood  perfusivity  and  tissue  conductivity  are  expected 
than  in  absorption  and  scattering  coefficients  in  the  cali- 
bration process. 

Using  nested-block  optimization  algorithm,  calibrations 
were  performed  using  different  sets  of  thermal  images  in 
the  objective  function  with  number  of  images  N = 10,  20, 
30,  40,  50,  60.  The  variation  in  the  calibrated  coefficients  is 
shown  in  Fig.  6a-d.  It  is  observed  that  the  perfusion 
coefficients  and  thermal  conductivity  parameters  are  the 
most  sensitive  to  the  calibration.  Figure  6e  and  f shows  the 
sensitivity  of  the  objective  function  and  its  gradient  to  the 
number  of  thermal  images  used  in  the  optimizations.  It  was 


observed  that  the  objective  function  is  rather  insensitive  to 
the  number  of  thermal  images  used  in  the  optimizations  of 
N > 30.  For  time  critical  applications,  this  implies  that  the 
calibration  process  may  be  just  as  accurate  using  only  half 
the  full  data  set.  For  comparison  at  a different  time 
instance,  thermal  contour  and  a cut-line  of  the  FEM 
predictions  using  N = 20  and  60  is  given  in  Fig.  7 at  time 
t = 235  s. 


5 Discussion  and  conclusions 

The  extent  of  thermal  damage  associated  with  traditional 
laser  therapies  is  limited  by  thermal  diffusion  in  the  tissue. 
Most  conventional  laser  therapies,  other  than  photody- 
namic therapy,  do  not  permit  selective  destruction  of  tumor 
tissue  while  preserving  the  surrounding  healthy  tissue.  In 
order  to  circumvent  the  limitations  of  traditional  laser 
therapies  employed  in  prostate  cancer  treatment,  nanoshells 
were  employed  to  enhance  the  thermal  deposition  and 
selectivity  of  thermal  damage. 

In  this  study,  a new  computational  methodology  for 
solving  the  Pennes  bioheat  transfer  equation  which  incor- 
porates a nested-block  optimization  algorithm  is  proposed. 
This  methodology  is  used  to  simulate  the  transient  tem- 
perature field  during  laser  surgery  on  a prostate  tumor 
inoculated  on  the  back  of  SCID  mice.  The  numerical 
simulations  for  this  test  case  are  in  good  agreement  with 
the  MRTI  data  from  the  surgical  procedure  but  required 
several  iterations  through  the  calibration  process.  For  the 
single  test  case  presented  herein  the  calibration  process 
stabilized  after  30  time  steps  (150  seconds  into  the 
surgery). 

A flexible  Monte  Carlo  based  formulation  for  the  flu- 
ence  in  the  laser  source  term  is  incorporated  in  the  Pennes 
bioheat  equation,  which  is  solved  by  a finite  element 
method.  This  formulation  is  considerably  more  general 
than  the  classical  approach  allowing  for  non-homoge- 
neous  and  anisotropic  tissues,  and  extends  the  range  of 
validity  to  include  more  complex  tumor  geometry.  Pre- 
liminary results  using  this  more  general  laser  model  show 
only  modest  variation  (lest  than  5%)  in  the  temperature 
field  when  compared  to  results  produced  using  a classical 
formulation  for  the  fluence.  We  anticipate  the  new 
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Fig.  5 Pre-calibration 
comparison  of  FEM  prediction 
to  thermal  images  at 
time  = 300  s.  a MRTI  data 
interpolated  onto  the  FEM 
mesh,  b unfiltered  MRTI  data 
and  filtered  MRTI  data,  c FEM 
prediction  using  isotropic  laser 
source  term  for  heating 
d filtered  MRTI  data  and  FEM 
prediction  used  isotropic  laser 
source  term  for  heating,  e FEM 
prediction  using  Monte  Carlo 
laser  source  term  for  heating, 
f filtered  MRTI  data  and  FEM 
prediction  used  Monte  Carlo 
laser  source  term  for  heating 


Table  1 Model  parameters  used  in  Pennes  bioheat  transfer  model 


us  14.2 

[cm”1] 

Ha  2.15  [cm  1 ] 

Ptissue 

1,045  [kg  m-3] 

Ta  [308  K] 

c“ood  3,840  [ J kg_1K_1 

i tissue 

J cp 

3,600 

[J  kg-'K-1] 

Table  2 Numerical  values  of  thermal  conductivity  and  blood  perfusivity 
used  in  Pennes  bioheat  transfer  model  at  the  beginning  of  calibration  [5] 


k(u) 

= kQ  + k\  atan (k2(u—k3)) 

0)(u)  = 

oj0  + co  i atan (co2(u—co3)) 

k0 

0.561  [J  s-’nr'tr1] 

CD0 

0.255  [k  g-V'nT1] 

k\ 

0.0427  [J  s^'nr'KT1] 

COi 

-0.137  [k  g-'s'Vr1] 

k2 

0.0252  [K-1] 

C02 

2.3589  [K_1] 

h 

315.0  [K] 

co3 

315.0  [K] 
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Fig.  6 Calibration  results  a-d  Calibrated  thermal  conductivity  and 
blood  perfusivity  parameters  as  a function  of  the  number  of  thermal 
images  used  in  the  objective  function,  N = 10,  20,  30,  40,  50,  60. 
Sensitivity  of  Objective  Function  e and  f The  value  of  the  objective 


(f)  e*io4 


20  22  24  26  28  30 

iteration  # 


function,  (5),  is  plotted  as  a function  of  the  number  of  thermal  images 
used  to  obtain  the  calibrated  parameters,  which  indicate  that  both 
objective  functional  and  its  gradient  become  insensitive  to  the  number 
of  images  included  in  the  calibration  if  N > 30 
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Fig.  7 Post-calibration 
comparison  of  FEM  prediction 
to  thermal  images  at 
time  = 235  s.  a MRTI  data 
interpolated  onto  the  FEM 
mesh,  b unfiltered  MRTI  data 
and  filtered  MRTI  data,  c 
Calibrated  FEM  prediction  with 
Monte  Carlo  laser  source  term 
for  heating  using  20  of  the  60 
available  thermal  images  for 
calibration,  d filtered  MRTI 
data  and  FEM  prediction  using 
Monte  Carlo  laser  source  term 
for  heating,  e Calibrated  FEM 
prediction  with  Monte  Carlo 
laser  source  term  for  heating 
using  50  of  the  60  available 
thermal  images  for  calibration, 
f filtered  MRTI  data  and  FEM 
prediction  using  Monte  Carlo 
laser  source  term  for  heating. 


(c)  310  I<  - 340  K 


(e)  310  K - 340  K 


formulation  will  be  a more  significant  factor  in  the 
analysis  as  we  investigate  more  complex  tumor  geometry 
from  the  perspective  of  laser  source  modeling.  Additional 
validation  tests  involving  more  specimens  with  different 


laser  treatments  protocols  are  underway  to  further  refine 
this  methodology. 

The  computer  model  developed  in  this  study  for  laser- 
induced  thermal  therapy  is  capable  of  predicting  the 
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temperature  due  to  laser  heating  in  a prostate  tumor  with 
nanoshell  inclusion  by  incorporating  changes  in  the 
absorption  and  scattering  effects.  Using  this  model,  optimal 
hyperthermia  protocols  can  be  developed.  With  maturity  of 
nano- technology  such  as  nanoshell  injection  and  a reliable 
computational  model,  the  outcome  of  laser  surgery  for 
cancer  treatment  will  become  more  controllable  and  pre- 
dictable. MRTI  provides  an  in  vivo  measurement  tool  for 
mapping  the  spatiotemporal  temperature  elevation  in  the 
tumor  during  and  following  laser  therapy.  This  technique 
for  temperature  measurement  has  resolution  of  less  than 
1°C.  With  the  nested-block  optimization  algorithm  pre- 
sented in  this  study,  tissue  properties  can  be  accurately 
characterized  during  the  course  of  a surgical  operation  and 
resulting  in  more  reliable  prediction  of  temperature  field 
and  treatment  outcomes. 

Acknowledgments  The  support  of  this  work  by  the  National  Sci- 
ence Foundation  under  grant  CNS-0540033  and  by  National  Institutes 
of  Health  under  Small  Animal  Imaging  Facility  core  grant  CA  16672 
are  gratefully  acknowledged. 


References 

1.  Baer  R,  Neuhauser  D,  Weiss  S (2004)  Enhanced  absorption 
induced  by  a metallic  nano  shell.  Nano  Fett  4(1):  85-88 

2.  Benson  SJ,  Mclnnes  FC,  More  J,  Sarich  J (2005)  TAO  user 
manual  (revision  1.8).  Technical  Report  ANF/MCS-TM-242, 
Mathematics  and  Computer  Science  Division,  Argonne  National 
Faboratory.  http://www.mcs.anl.gov/tao 

3.  Brigger  I,  Dubernet  C,  Couvreur  P (2002)  Nanoparticles  in  cancer 
therapy  and  diagnosis.  Adv  Drug  Deliv  Rev  54(5) :63 1-51 

4.  Demkowicz  F,  Rachowicz  W,  Devloo  Ph  (2002)  A fully  auto- 
matic hp- adaptivity.  J Sci  Comput  17(1-3):  127-155 

5.  Diller  KR,  Valvano  JW,  Pearce  JA  (2005)  Bioheat  transfer.  In: 
Kreith  F,  Goswami  Y (eds)  The  CRC  handbook  of  mechanical 
engineering,  2nd  edn.  CRC  Press,  Boca  Raton,  pp  4-278-4-357 

6.  Feng  Y,  Rylander  MN,  Bass  J,  Oden  JT,  Diller  K (2005)  Optimal 
design  of  laser  surgery  for  cancer  treatment  through  nanoparticle- 
mediated  hyperthermia  therapy.  In:  NSTI-Nanotech  2005,  vol  1. 
pp  39-42 

7.  Ferrari  M (2005)  Cancer  nanotechnology:  opportunities  and 
challenges.  Nat  Rev  Cancer  5(3):  161-171 

8.  Hindman  JC  (1966)  Proton  resonance  shift  of  water  in  the  gas  and 
liquid  states.  J Chem  Phys  44:4582-4592 

9.  Hirsch  FR,  Stafford  RJ,  Bankson  JA,  Sershen  SR,  Rivera  B,  Price 
RE,  Hazle  JD,  Halas  NJ,  West  JF  (2003)  Nanoshell-mediated 
near-infrared  thermal  therapy  of  tumors  under  magnetic  reso- 
nance guidance.  Proc  Natl  Acad  Sci  100(23):  13549-13554 


10.  Jain  PK,  El-Sayed  IH,  E-Sayed  MA  (2007)  Au  nanoparticles 
target  cancer,  nanotoday 

11.  Fin  A,  Hirsch  F,  Fee  MH,  Barton  J,  Halas  N,  West  J,  Drezek  R 
(2004)  Nanoshell-enabled  photonics-based  imaging  and  therapy 
of  cancer.  Technol  Cancer  Res  Treat  3(1) 

12.  Masters  A,  Steger  AC,  Fees  WR,  Walmsley  KM,  Bown  SG 
(1992)  Interstitial  laser  hyperthermia:  a new  approach  for  treating 
liver  metastases.  Br  J Cancer  66(3) :5 18-522 

13.  Nagata  Y,  Hiraoka  M,  Akuta  K,  Abe  M,  Takahashi  M,  Jo  S, 
Nishimura  Y,  Masunaga  S,  Fukuda  M,  Imura  H (1990)  Radio- 
frequency thermotherapy  for  malignant  liver  tumors.  Cancer 
65(8):  1730-1736 

14.  Nau  WH,  Roselli  RJ,  Milam  DF  (1999)  Measurement  of  thermal 
effects  on  the  optical  properties  of  prostate  tissue  at  wavelengths 
of  1,  064  and  633  nm.  Fasers  Surg  Med  24(l):38-47 

15.  Niemz  MH  (2004)  Faser- tissue  interactions:  fundamentals  and 
applications,  3rd  edn.  Springer,  Berlin 

16.  Oden  JT,  Diller  KR,  Bajaj  C,  Browne  JC,  Hazle  J,  Babuska  I, 
Bass  J,  Demkowicz  F,  Feng  Y,  Fuentes  D,  Prudhomme  S, 
Rylander  MN,  Stafford  RJ,  Zhang  Y (2007)  Dynamic  data-driven 
finite  element  models  for  laser  treatment  of  prostate  cancer.  Num 
Meth  PDE  23:904-922 

17.  Olsrud  J,  Wirestam  R,  Brockstedt  S,  Nilsson  AMK,  Tranberg 
KG,  Stahlberg  F,  Persson  BRR  (1998)  MRI  thermometry  in 
phantoms  by  use  of  the  proton  resonance  frequency  shift  method: 
application  to  interstitial  laser  thermotherapy.  Phys  Med  Biol 
43(9):2597-2613 

18.  Pennes  HH  (1948)  Analysis  of  tissue  and  arterial  blood  temper- 
atures in  the  resting  forearm.  J Appl  Physiol  1:93-122 

19.  Robinson  DS,  Parel  JM,  Denham  DB,  Gonzalez- Cirre  X,  Manns 
F,  Milne  PJ,  Schachner  RD,  Herron  AJ,  Comander  J,  Hauptmann 
G (1998)  Interstitial  laser  hyperthermia  model  development  for 
minimally  invasive  therapy  of  breast  carcinoma.  J Am  Coll  Surg 
186(3):284-292 

20.  Rylander  MN,  Feng  Y,  Bass  J,  Diller  KR  (2005)  Thermally 
induced  injury  and  heat-shock  protein  expression  in  cells  and 
tissues.  Ann  N Y Acad  Sci  1066:222-242 

21.  Rylander  MN,  Feng  Y,  Zhang  J,  Bass  J,  Stafford  RJ,  Hazle  J, 
Diller  K (2006)  Optimizing  hsp  expression  in  prostate  cancer 
laser  therapy  through  predictive  computational  models.  J Biomed 
Opt  11:4:041113 

22.  U.S.  National  Institutes  of  Health.  National  Cancer  Institute. 
http://www.cancer.gov 

23.  Welch  AJ,  van  Gemert  MJC  (1995)  Optical-thermal  response  of 
laser-irradiated  tissue.  Plenum  Press,  New  York 

24.  Welsh  AJ,  Gardner  C (2002)  Optical  and  thermal  response  of 
tissue  to  laser  radiation.  In:  Waynant  RW  (ed)  Fasers  in  medi- 
cine, Chap.  2.  pp  27-45 

25.  Wust  P,  Hildebrandt  B,  Sreenivasa  G,  Rau  B,  Gellermann  J,  Riess 
H,  Felix  R,  Schlag  PM  (2002)  Hyperthermia  in  combined  treat- 
ment of  cancer.  Fancet  Oncol  3(8):487-497 

26.  Yang  R,  Sanghvi  NT,  Rescorla  FJ,  Kopecky  KK,  Grosfeld  JF 
(1993)  Fiver  cancer  ablation  with  extracorporeal  high-intensity 
focused  ultrasound.  Eur  Urol  23(1):  17-22 


l 


Springer 


frecpopernne 


